clear
clc

parfor zz=1:68
tic    
%%%%%loading data    
w=load('competitions.mat');
COMP=w.comp_list(zz,1);
x=load('data.mat');
y=load('qestimates.mat');
competitionid=x.competitionid;
maxscore = x.maxscore(competitionid==COMP,1);
pubscore_normal = x.pubscore_normal(competitionid==COMP,1);
length_days= x.length_days(competitionid==COMP,1);
length_days=length_days(1);
rewardquantity= x.rewardquantity(competitionid==COMP,1);
rewardquantity=rewardquantity(1);
mergers_data= x.totmergers(competitionid==COMP,1);
mergers_data=mergers_data(1);
nSim = 500;

%%%%%state space parameters
N=40; %number of individuals
T=360; %number of periods
Svalues=unique([unique(maxscore);max(maxscore)+[0.002:0.002:0.08]']);
nScore = 20;
Svalues = Svalues(3:nScore+2,1);
nsubs_data=size(maxscore,1);
lambda1=0.75;
lambda2=1-lambda1;

%%%%%q-function
betachangehat=[y.beta0(y.cid==COMP,1) ; y.beta1(y.cid==COMP,1); y.beta1tr(y.cid==COMP,1)];
probchange = @(score) (exp(betachangehat(1) + betachangehat(2)*score)./(1+exp(betachangehat(1) + betachangehat(2)*score)));
q_sp=probchange(Svalues);
%q_sp(end,1)=0;
probchange = @(score) (exp(betachangehat(1) + betachangehat(3) + betachangehat(2)*score)./(1+exp(betachangehat(1)+ betachangehat(3) + betachangehat(2)*score)));
q_team=probchange(Svalues);
%q_team(end,1)=0;

%%%%%estimation

theta0=[0.01,0.4;0.01,0.6; 0.08,0.4;0.08,0.6;([log(nsubs_data/(T*lambda1))/log(0.0025),log(mergers_data/(T*lambda2))/log(0.005)])];
f0 = zeros(size(theta0,1),1);
for r=1:size(theta0,1)
   [mergers,subs,maxs,~]=equilibrium(theta0(r,1),theta0(r,2),N,T,lambda1,lambda2,q_team,q_sp,nScore,Svalues,mergers_data,nSim);
    f0(r,1)=((mergers-mergers_data)/mergers_data)^2 + ((subs-nsubs_data)/nsubs_data)^2;
end
[~ , i] = min(f0);
theta0=log(theta0(i,:));

%optimization
options=optimset('MaxIter',500,'MaxFunEval',5000,'Display','iter','TolX',1e-6,'TolFun',1e-4);
[thetahat, fhat] = fminsearch(@(x) objfun(x,N,T,lambda1,lambda2,q_team,q_sp,nScore,Svalues,mergers_data,nsubs_data,nSim),theta0,options);

sigmaS = exp(thetahat(1));
sigmaT = exp(thetahat(2));
zz
fhat
[sigmaS,sigmaT]

%%%%%fit and counterfactuals
[mergers,subs,maxs,Omega]=equilibrium(sigmaS,sigmaT,N,T,lambda1,lambda2,q_team,q_sp,nScore,Svalues,mergers_data,nSim);
fit = [mergers_data,mergers,nsubs_data, subs,max(maxscore), maxs];
fit
 
% 
% options_g=optimset('MaxIter',1,'MaxFunEval',10,'Display','off');    
% [~,~,~,~,GRAD,~] = fminunc(@(x) objfun(x,N,T,lambda1,lambda2,q_team,q_sp,nScore,Svalues,mergers_data,nsubs_data,nSim),thetahat,options_g) 

[mergers_nm,subs_nm,maxs_nm]=equilibrium(sigmaS,1e+5,N,T,lambda1,lambda2,q_team,q_sp,nScore,Svalues,mergers_data,nSim);
counter=[mergers,mergers_nm,subs,subs_nm,maxs,maxs_nm];
% counter

parsave(sprintf('Estimates/%s_%02d.mat', 'estimates',COMP),sigmaS,sigmaT,COMP,N,T,fit,counter,fhat,lambda1,lambda2,q_team,q_sp,nScore,Svalues,mergers_data,nsubs_data,nSim)
toc
end